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Motivated by the need to statistically quantify the difference between two spatio-temporal 
datasets that arise in climate downscaling studies, we propose new tests to detect the differences 
of the covariance operators and their associated characteristics of two functional time series. Our 
two sample tests are constructed on the basis of functional principal component analysis and 
self-normalization, the latter of which is a new studentization technique recently developed for 
the inference of a univariate time series. Compared to the existing tests, our SN-based tests 
allow for weak dependence within each sample and it is robust to the dependence between 
the two samples in the case of equal sample sizes. Asymptotic properties of the SN-based test 
statistics are derived under both the null and local alternatives. Through extensive simulations, 
our SN-based tests are shown to outperform existing alternatives in size and their powers are 
found to be respectable. The tests are then applied to the gridded climate model outputs and 
interpolated observations to detect the difference in their spatial dynamics. 

Keywords: climate downscaling; functional data analysis; long run variance matrix; 
self-normalization; time series; two sample problem 


1. Introduction 

Functional data analysis (FDA) which deals with the analysis of curves and surfaces has 
received considerable attention in the statistical literature during the last decade (Ram¬ 
say and Silverman [19, 20] and Ferraty and Vieu [6]). This paper falls into a sub-field 
of functional data analysis: inference for temporally dependent functional data. Specifi¬ 
cally, we focus on testing the equality of the second-order structures (e.g., the covariance 
operators and their associated eigenvalues and eigenfunctions) of two temporally depen¬ 
dent functional sequences. Our work is partially motivated by our ongoing collaboration 
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with atmospheric scientists on the development and assessment of high-resolution cli¬ 
mate projections through statistical downscaling. Climate change is one of the most 
urgent problems facing the world this century. To study climate change, scientists have 
relied primarily on climate projections from global/regional climate models, which are 
numerical models that involve systems of differential equations and produce outputs at 
a prespecified grid. As numerical model outputs are widely used in situations where real 
observations are not available, it is an important but still open question whether the nu¬ 
merical model outputs are able to mimic/capture the spatial and temporal dynamics of 
the real observations. To partly answer this question, we view the spatio-temporal model 
outputs and real observations as realizations from two temporally dependent functional 
time series defined on the two-dimensional space and test the equality of their second- 
order structures which reflects their spatial dynamics/dependence. 

Two sample inference for functional data has been investigated by a few researchers. 
Fan and Lin [5], Cuevas et al. [4] and Horvath et al. [10] developed the tests for the equal¬ 
ity of mean functions. Benko et al. [1], Panaretos et al. [16], Fremdt et al. [7], and Kraus 
and Panaretos [12] proposed tests for the equality of the second-order structures. All 
the above-mentioned works assumed the independence between the two samples and/or 
independence within each sample. However, the assumption of independence within the 
sample is often too strong to be realistic in many applications, especially if data are 
collected sequentially over time. For example, the independence assumption is question¬ 
able for the climate projection data considered in this paper, as the model outputs and 
real station observations are simulated or collected over time and temporal dependence 
is expected. Furthermore the dependence between numerical model outputs and station 
observations is likely because the numerical models are designed to mimic the dynamics 
of real observations. See Section 5 for empirical evidence of their dependence. In this 
paper, we develop new tests that are able to accommodate weak dependence between 
and within two samples. Our tests are constructed on the basis of functional principal 
component analysis (FPCA) and the recently developed self-normalization (SN) method 
(Shao [21]), the latter of which is a new studentization technique for the inference of a 
univariate time series. 

FPCA attempts to find the dominant modes of variation around an overall trend 
function and has been proved a key technique in the context of FDA. The use of FPCA 
in the inference of temporally dependent functional data can be found in Gabrys and 
Kokoszka [8], Hormann and Kokoszka [9], Horvath et al. [10] among others. To account 
for the dependence, the existing inference procedure requires a consistent estimator of 
the long run variance (LRV) matrix of the principal component scores or consistent 
estimator of the LRV operator. However, there is a bandwidth parameter involved in 
the LRV estimation and its selection has not been addressed in the functional setting. 
The same issue appears when one considers the block bootstrapping and subsampling 
schemes (Lahiri [13] and Politis et al. [18]), since these techniques also require the selection 
of a smoothing parameter, such as the block length in the moving block bootstrap, 
and the window width in the subsampling method (see, e.g., Politis and Romano [17] 
and McMurry and Politis [15]). Since the finite sample performance can be sensitive 
to the choice of these tuning parameters and the bandwidth choice can involve certain 
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degree of arbitrariness, it is desirable to use inference methods that are free of bandwidth 
parameters. To this end, we build on the bandwidth-free SN method (Shao [21]) recently 
developed in the univariate time series setup, and propose SN-based tests in the functional 
setting by using recursive estimates obtained from functional data samples. 

In time series analysis, the inference of a parameter using normal approximation typ¬ 
ically requires consistent estimation of its asymptotic variance. The main difficulty with 
this approach (and other block-based resampling methods) is the sensitivity of the finite 
sample performance with respect to the bandwidth parameter, which is often difficult 
to choose in practice without any parametric assumptions. As a useful alternative, the 
self-normalized approach uses an inconsistent bandwidth-free estimate of asymptotic vari¬ 
ance, which is proportional to asymptotic variance, so the studentized quantity (statistic) 
is asymptotically pivotal. Extending the early idea of Lobato [14], Shao [21] proposed a 
very general kind of self-normalizers that are functions of recursive estimates and showed 
the theoretical validity for a wide class of parameters of interest. The settings in the 
latter two papers are however limited to univariate time series. The generalization of the 
SN method from univariate to functional time series was first done in Zhang et al. [23] 
where the focus was on testing the structure stability of temporally dependent functional 
data. Here we extend the SN method to test the equality of the second-order properties 
of two functional time series, which is rather different and new techniques and results are 
needed. To study the asymptotic properties of the proposed test statistics, we establish 
functional central limit theorems for the recursive estimates of quantities associated with 
the second-order properties of the functional time series which seems unexplored in the 
literature and are thus of independent interest. Based on the functional central limit 
theorem, we show that the SN-based test statistics have pivotal limiting distributions 
under the null and are consistent under the local alternatives. From a methodological 
viewpoint, this seems to be the first time that the SN method is extended to the two 
sample problem. Compared to most of the existing methods which assumed the indepen¬ 
dence between the two samples and/or independence within each sample, the SN method 
not only allows for unknown dependence within each sample but also allows for unknown 
dependence between the two samples when the sample sizes of the two sequences are 
equal. 


2. Methodology 

We shall consider temporally dependent functional processes {(Aj(t), Yi(t)), t £ X}ff?f 
defined on some compact set X of the Euclidian space, where X can be one-dimensional 
(e.g., a curve) or multidimensional (e.g., a surface or manifold). For simplicity, we consider 
the Hilbert space Et of square integrable functions with I = [0,1] (or X = [0, l] 2 ). For any 
functions f,g £ EL the inner product between / and g is defined as J x f(t)g(t) dt and 
|| • || denotes the inner product induced norm. Assume the random elements all come 
from the same probability space (f2,A, V). Let L p be the space of real valued random 
variables with finite L p norm, that is, (E |A| p ) 1 / p < oo. Further, we denote L ^ the space 
of H valued random variables X such that (£l||A’|| p ) 1 / p < oo. 


4 


X. Zhang and X. Shao 


Given two sequences of temporally dependent functional observations, {Xi{t)}^l 1 and 
{Yi(t)}^ 1 defined on a common region X, we are interested in comparing their second- 
order properties. Suppose that the functional time series are second-order stationary. We 
assume that E[Xi(t)] = E[Yi(t)] = 0. The result can be easily extended to the situation 
with nonzero mean functions. Define Cx = E[{Xi, -)Xj\ and Cy = E[(Yi, •)!*] as the co- 
variance operators of the two sequences respectively. For the convenience of presentation, 
we shall use the same notation for the covariance operator and the associated covariance 
function. Denote by {</>x}j2=i an d the eigenfunctions and eigenvalues of Cx- 

Analogous quantities are and {Ay}^ for the second sample. Denote by |v| the 

Euclidean norm of a vector v £ R p . Let vech(-) be the operator that stacks the columns 
below the diagonal of a symmetric m x m matrix as a vector with m{m + l )/2 com¬ 
ponents. Let £>[0,1] be the space of functions on [0,1] which are right-continuous and 
have left limits, endowed with the Skorokhod topology (see Billingsley [2]). Weak conver¬ 
gence in D[0,1] or more generally in the R m -valued function space D m [0,l] is denoted 
by “=>”, where m £ N and convergence in distribution is denoted by > d ”. Define [aj 
the integer part of a £ R, and Sij = 1 if i = j and Sij =0 if i ft j ■ In what follows, we 
shall discuss the tests for comparing the three quantities Cx, ftx and ^x w ith Cy, ft Y 
and Ay, respectively. 

2.1. Covariance operator 

Consider the problem of testing the hypothesis H\q : Cx = Cy versus the alternative 
H[ a ■ Cx ft Cy (in the operator norm sense) for two mean zero stationary functional 
time series {X^t)}^ and Let N = Ni + N 2 . Throughout the paper, we 

assume that 


fVi/X-^ 71 , N 2 /N ^ 72 , as min(IVi, N 2 ) ->■ +oo, 

where 71,72 £ (0,1) and 71 +72 = 1. Define the one-dimensional operator A, = ( Xi , -)Xi = 
Xi ® Xi and y 3 ={Yj,-)Yj = Y 3 ®Yj. Let Cxy be the empirical covariance operator based 
on the pooled samples, that is, 


/ n 1 n 2 \ 

< 2I > 

Denote by {A^y } and {4>xy} the corresponding eigenvalues and eigenfunctions. The pop¬ 
ulation counterpart of Cxy is then given by Cxy — 71 Cx +72 Cy whose eigenvalues and 
eigenfunctions are denoted by {A J } and {ft} respectively. Further let Cx,m = 77 YftiLi %% 
be the sample covariance operator based on the subsample {Xi(t)}^L 1 with 2 < m < N±. 
Define {ft J x m }JLi and {Ay m }JLi the eigenfunctions and eigenvalues of Cx,m respec¬ 
tively, that is, 

J^Cx,m(tc)ftx, m ( s ) ds = %, m ^f >m (t), 


( 2 . 2 ) 
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and J I $ic tm {t)^ m (t)dt = 5i j . Similarly, quantities C Y , m >, {ft YtTn ,}fl i and {Ay^/}^ 
are defined for the second sample with 2 < m' < N 2 - To introduce the SN-based test, we 
define the recursive estimates 


c fe J — {{Cx, [fciVi /n\ — Cy, [kN 2 /n\ )<t> l XYi$x Y )i 2 < k < N, 1 < i , j < K , 

which estimate the difference of the covariance operators on the space spanned by 
■ Here K is a user-chosen number, which is held fixed in the asymptotics. Denote 

by &k = vech(Cfe) with C= (c^, J )f J=1 . In the independent and Gaussian case, Panaretos 
et al. [16] proposed the following test (hereafter, the PKM test), 


Tn 1} n 2 


2N hh ’ 




W XY! 




m^xv)f 


which converges to x(k+i)k /2 un der the null. To take the dependence into account, we 
introduce the SN matrix 

1 N 

Vsn,n( d) = & (&k ~ div)(dfc — ax) , (2-3) 

fc= 1 

with d=(K + l)K/2. The SN-based test statistic is then defined as, 

G ( il tN (d) = Na'xiV^id))- 1 ^. (2.4) 

Notice that the PKM test statistic can also be written as a quadratic form of djv but 
with a different normalization matrix that is only applicable to the independent and 
Gaussian case. The special form of the SN-based test statistic makes it robust to the 
dependence within each sample and also the dependence between the two samples when 
their sample sizes are equal. We shall study the asymptotic behavior of Gg 1 ^- N (d) under 
the weak dependence assumption in Section 3. 


2.2. Eigenvalues and eigenfunctions 

In practice, it is also interesting to infer how far the marginal distributions of two se¬ 
quences of stationary functional time series coincide/differ and quantify the difference. 
By the Karhunen-Loeve expansion (Bosq [3], page 26), we have 



i=i j= i 


where Px t ,j = f x Xi(t)</>x(t) d t and j3 Yi ,j = f x li(f)<^y(t) d t are the principal components 
(scores), which satisfy that E\^Xi,jPxi,j'] = Sjf and E[/3 Yi ,jf3 Yi ,j'] =Sjj>- The prob¬ 
lem is then translated into testing the equality of the functional principal components 
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(FPC’s) namely the eigenvalues and eigenfunctions. For a prespecified positive inte¬ 
ger M, we denote the vector of the first M eigenvalues by A^ M = (A^-,..., A]Sf) and 
Xy M = (Ay,..., Xy)- Further define 4>x M = (4>x,- ■ ., <f>x) and </>y M = (4>y, ■ ■ ■ > (fry ) the 
first M eigenfunctions of the covariance operators Cx and Gy, respectively. Since the 
eigenfunctions are determined up to a sign, we assume that (0^-. 4>y) > 0 in order for the 
comparison to be meaningful. We aim to test the null hypothesis # 2,0 : ^x M = A y M and 
# 3,0 : </>; 1/1 = </>y M versus the alternatives that H 2t a ■ A \ M ^ Ay M and H 3ja : ^ </>y M 

(in the L 2 norm sense). The problem of comparing the FPC’s of two independent and 
identically distributed (i.i.d.) functional sequences has been considered in Benko et al. 
[1], where the authors proposed an i.i.d. bootstrap method which seems not applicable 
to the dependent case. The block bootstrap based method is expected to be valid in the 
weakly dependent case but the choice of the block size seems to be a difficult task in the 
current setting. To accommodate the dependence and avoid the bandwidth choice, we 
adopt the SN idea. 

Recall the recursive estimates of the eigenvalues X J x m and Ay m , which are calculated 
based on the subsamples {X^t)}^ and {^(t)}^. Let 0 J k = X 3 x>[kNl/Ni ~ ^ Y ,[kN 2 /N\ 

and 9 k = 0\,Off)' with [Ne J <k<N for some e € (0,1], which is held fixed in the 
asymptotics. We consider the trimmed SN-based test statistic 


( N 

g { sn, n ( m ) = n3 0n\ k 2 0 k -6 N )0 k -d N y 

lfc=[iVeJ 



(2.5) 


The trimmed version of the SN-based test statistic is proposed out of technical consider¬ 
ation when the functional observations lie on an infinite dimensional space. It can be seen 
from the proof in the supplemental material [22] that the trimming is not required when 
functional data lie on a finite-dimensional space; see Remark 0.1 in the supplemental 
material [22]. 


Remark 2.1. To compare the difference between the eigenvalues, one may also con¬ 
sider their ratios. Define ( k = (^x : i k N 1 /N\/^y,[kN 2 /N\^ ■ ■ ■ i^x,lkN 1 /N\/^Y,[kN 2 /N\y f° r 
k = L-NeJAn alternative SN-based test statistic is given by 

— 1 m)'(^ fc 2 (c*- 6 v)(a- 6 v)'} ( 6 v- 1 m), (2.6) 

l fc=[jVeJ J 


where 1 m is a M-dimensional vector of all ones. Since the finite sample improvement 
by using Gg^ N (M) is not apparent, we do not further investigate the properties of 

Gsn,n(M)- 

We now turn to the problem of testing the equality of the eigenfunctions. To proceed, 
we let 


Vj = 0Xv > fix y >■■■■> ftxy) 


( 2 . 7 ) 
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be a vector of p—j orthonormal basis functions for j = 1 , 2 ,..., M with M <p and p being 
a user chosen number. Recall that (p 3 x m (t) and (pf. m , (t) are the jth eigenfunctions of the 

empirical covariance operators Cx,m and CV,m' which are computed based on the first 
m (and to') samples. Here we require that (<f? x m ,(f> x N } > 0 and (ftym'trx n x ) — 0 f° r 
2 <m < Ni and 2<m!< N% ■ As the eigenfunctions are defined on an infinite-dimensional 
space, we project the difference between the jth eigenfunctions onto the space spanned 
by Uj. Formally, we define the projection vectors 

Vk = ((&x,\.kNi/N] ~ &Y,[kN 2 /N\ > fe)’ • ■ • > (^X,|_fciVi/iVJ “ ^Y,\kN 2 /N\ > fc))> 

where 1 < j < M and k — |_-/VeJFurther let fjk = (.Vk’Vk’ • • • >£ ® M ° with 
M 0 = M ( 2 P~ j ' / ~ 1 ) _ The trimmed SN-based test statistic is then defined as 

g< sn,n( m o)= n VnIjp ^2 k 2 (f) k -fj N )(f) k -rj N )'\ fj N , ( 2 . 8 ) 

l fc=LJVeJ J 

for some 0 < e < 1. 


Remark 2.2. It is worth noting that N (M, o) is designed for testing the equality 

of the first M eigenfunctions. Suppose we are interested in testing the hypothesis for a 
particular eigenfunction, that is, the null <f? x = (jyL versus the alternative (f) J x ^ (f> Y . We 
can consider the basis functions 


Fj — ($xy ■ ■ ■ > &xy i &xy ^xy)’ 


and the projection vector rf k = (((f> x kN . t /n \ ' 

^Y,[kN 2 /N\ ’ yXy)i (^X,LfcJVi/JVJ ~ iY,[kN 2 /N\ ’ YXy) • ‘ • ’ ( < ^X,\kNi/N\ ~ VY,\kN 2 /N\ > VXY/ 

The SN-based test statistic can then be constructed in a similar manner. We also note 


$Y, \_kN 2 /N\ > $xy)i ■i (fix, [fcJVi/JVJ 

A P xy)Y- 


that when <jP x ^ cf> Y and <p x = cj) Y for i ^ j, the choice of Vj may result in trivial power 
because (<j? x — </>y,</> 2 ) for i ^ j can be close to 0. In this case, one remedy is to consider 
alternative basis functions, for example, (4.5) and (4.6) as suggested in the simulation. 


Remark 2.3. The choice of the basis functions Oj is motivated by the Bahadur rep¬ 
resentation of the recursive estimates in the supplemental material [22]. Under suitable 
assumptions as given in the next section, it can be shown that 


! </*) — (Y’x > 4>) ~ ^ 

i=1 


s^a 


I3xi,s/3xi,c 

»x - K 


■ i </*) r + 


:,ki 


(2.9) 


with R x k being the remainder term and 4> £ Lr{T j. The second term on the RHS of (2.9) 
plays a key role in determining the limiting distribution of the SN-based test statistic. 
When = with j ^ a, the linear term reduces to — t Y2l=i ■ which satisfies 

the functional central limit theorem under suitable weak dependence assumption. Notice 
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that the linear term vanishes when = <p° x and the asymptotic distribution of the pro¬ 
jection vector is degenerate. It is also worth noting that the linear terms in the Bahadur 
representations of (fxk^° x ) an d Wxk’^x) are opposite of each other which suggests 
that when testing the eigenfunctions jointly, the basis functions should be chosen in a 
proper way so that the asymptotic covariance matrix of the projection vector, that is, fjk 
is nondegenerate. 


3. Theoretical results 

To study the asymptotic properties of the proposed statistics, we adopt the dependence 
measure proposed in Hormann and Kokoszka [9], which is applicable to the temporally 
dependent functional process. There are also other weak dependence measures (e.g., 
mixing) or specific processes (e.g., functional linear processes) suitable for the asymp¬ 
totic analysis of functional time series (see Bosq [3]), we decide to use Hormann and 
Kokoszka’s i p -m-approximating dependence measure for its broad applicability to linear 
and nonlinear functional processes as well as its theoretical convenience and elegance. 

Definition 3.1. Assume that {Xi\ £ L ^ with p > 0 admits the following representation 

Xi = f(£i,£i- 1 ,...), * = 1 , 2 ,..., (3.1) 

where the £i’s are i.i.d. elements taking values in a measurable space S and f is a mea¬ 
surable function f: S°° — > H. For each i £ N, let {ejd} JgZ be an independent copy of 
{£j}j&- The sequence {Xi} is said to be L p -m-approximable if 

OO 

J2mX m -X^\\ p ) 1/P <oo, (3.2) 

771=1 

where x\ m) = /(£*,£*_!,..• • ■)■ 

Define B q (r) as a g-dimensional vector of independent Brownian motions. For e £ [0,1), 
we let 

W q {e) = B q (l)'J q (e)~ 1 B q (l), where J q (e) = J (B q (r)-rB q (l))(B q (r)-rB q (l))'dr. 

The critical values of W q := W q (0) have been tabulated by Lobato [14]. In general, the 
quantiles of W q (e) can be obtained via simulation. To derive the asymptotic properties 
of the proposed tests, we make the following assumptions. 

Assumption 3.1. Assume ^ anc ^ 'iTWIitu — are both L 4 -m- 

approximable and they are mutually independent. 

Assumption 3.2. Assume {{Xi(t),Yi(t))}fX^ Clg xl is an L 4 -m-approximable se¬ 
quence. 
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Assumption 3.3. Assume Ay > X' 2 x > • • • > Ay 0+1 and Ay > Ay > • • • > Ay 0+1 , for 
some positive integer too > 2. 

Note that Assumption 3.2 allows dependence between {W(t)} and (Y;(f)}, which is 
weaker than Assumption 3.1. To investigate the asymptotic properties of G^ x N (d) un " 
der the local alternatives, we consider the local alternative H\ a :Cx — Cy = LC/y/N 
with C being a Hilbert-Schmidt operator, where L is a nonzero constant. Define 
A = ((Cy>*,^ J ))f^- =1 £ R KxK as the projection of C onto the space spanned by 
..., <f> K } and assume that vech(A) ^ 0 £ R d . The following theorem states the 
asymptotic behaviors of G^l, N (d) under the null and the local alternatives. 

Theorem 3.1. Suppose Assumptions 3.1, 3.3 hold with mo > K. Further assume that the 
asymptotic covariance matrices A* d (A* d )' given in Lemma 0.3 is positive definite. Then 
under H lfi , G^ N (d) ~^ d W d and under Hy a , lim\ L \_> +00 lim N ^. +00 G^^d) = +oo. 
Furthermore, if ji = 72 , then the conclusion also holds with Assumption 3.1 replaced by 
Assumption 3.2. 

It is seen from Theorem 3.1 that G^y x(d) has pivotal limiting distributions under the 
null and they are consistent under the local alternatives as L — > + 00 . It is worth noting 
that in our asymptotic framework, d (or K) is assumed to be fixed as n — > 00 . Since K is 
usually chosen to make the first K principle components explain a certain percentage of 
variation (say 85%) , the magnitude of K critically depends on the prespecified threshold 
and the decay rate of the eigenvalues. In some cases, d = (K + \)K/2 can be quite 
large relative to sample size so it may be more meaningful to use the asymptotic results 
established under the framework that d —> 00 but d/n —> 0 as n —> 00 . This motivates the 
question that whether the following convergence result 

sup|P(Gg 1 ^ r N (d) < x) — P(Wd < x)| —> 0 as n —> 00 

iei 

holds when d diverges to 00 but at a slower rate than n. This would be an interesting 
future research topic but is beyond the scope of this paper. 

To study the asymptotics of Ggf, N {M) and Ggfj N (M 0 ), we introduce some notation. 
Let Wx. = /3xi,jPxi,k and k (h) = £’[(w^. — SjkXj)(uj : ’ x k +h — Sj'k'^j')] be the cross¬ 
covariance function between and at lag h. Set r x (h) := r J x ^ k (h). Define v x . = 
ui x . — E[uj x ] =co x — 8jk\j. Analogous quantities Ty ’ J ' k ( h) and Vy. can be defined for 
the second sample. We make the following assumption to facilitate our derivation. 

Assumption 3.f. Suppose that 

( +00 \ 2 +00 

J2 \ r x’ J ' k \ h )\) <+ 00 . J2 J2 \ r x( h )\<+°° 

h=—o o / j,k h=—oo 


(3.3) 
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and 


EE E I cum Kl> ’ v x H ’ v Xi 2 > v Xi 3 ) I < 00 • 

j,k j',k' *1 ,i 2 ,i 3 £Z 

TTie summability conditions also hold for the second sample {Yi(t)}. 


(3.4) 


Assumption 3.4 is parallel to the summability condition considered in Benko et al. [1] 
(see Assumption 1 therein) for i.i.d. functional data. It is not hard to verify the above 
assumption for Gaussian linear functional process (see, e.g., Bosq [3]), as demonstrated 
in the following proposition. 

Proposition 3.1. Consider the linear process X,(t ) = bj£i—j(t), where £j(t) = 
with {zij} being a sequence of independent standard normal random 
variables across both index i and j. Let ir(h) = JA bibi+h- Assume that Aj < 00 ani ^ 

12h l 7r (^-)l < 00 ■ Then Assumption 3.f holds for {Ai(t)}. 

Theorem 3.2. Suppose Assumptions 3.1, 3.3, 3.f hold with mo > M and the asymptotic 
covariance matrix AmA' m given in Lemma 0.5 is positive definite. Then under # 2,0 > we 
have G^ x n (M) —> d Wm(U)- Under the local alternative a '■ A X M — X Y M = -4= A with 
A ^ 0 € R M , we have lim^i^^ limjv^+oo G^ NN (M) = + 00 . 


(ON 

In order to study the asymptotic properties of G y sx N (M 0 ) under the null and local 


alternative, we further make the following assumption. 

Assumption 3.5. Let = f xj m \t)<f J x (t) dt, where x[ m ' ) is the m-dependent ap¬ 
proximation of Xi(t) (see Definition 3.1). Suppose one of the following conditions holds: 


OO OO 


EE{ £ (fc»i-fi) 4 } 1/4 <“. 


m—1 j =1 


1=1 


1/4 


< 00 , 


(3.5) 


+OO 


Ek^,< 


< + 00 , 


2 <j<p. 


s= 1 


The same condition holds for the second sample {Y)(t)}. 


(3.6) 


Theorem 3.3. Suppose Assumptions 3.1, 3.3, 3.f and 3.5 hold with mo > M and the 
asymptotic covariance matrix Am 0 A' m given in Lemma 0.7 is positive definite. Then 

under H 3j0 , we have G ( g XN (M 0 ) —t d WM 0 (e). 

Proposition 3.2. Define A by replacing 4> x Xl , (f> Y N2 and (jr XY with t)P x , <f? Y and fA 
in the definition of t)n ■ Consider the local alternative H 3 a : A = Lxf/'/N with ff> ^ 0 £ 
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R M ° . Suppose Assumptions 3.1, 3.3, 3-4 and 3.5 hold with mo > M and the asymptotic 
covariance matrix Am 0 A' Mq given in Lemma 0.7 is positive definite. Then we have 

m™ A / 1 ”! 1 G SN,n( M o)=+°o 

\L\—>oo iv—»+oo ’ 


under Ho a . 

It is worth noting that the conclusions in Theorem 3.2, Theorem 3.3 and Proposition 3.2 
also hold with Assumption 3.1 replaced by Assumption 3.2 and 71 = 72 - Finally, we point 
out that condition (3.5) can be verified for Gaussian linear functional process as shown 
in the following proposition. 

Proposition 3.3. Consider the Gaussian linear process in Proposition 3.1. Assume that 
E,“i «"<< < 00. Then Assumption 3.4 and condition (3.5) 

are satisfied for {Xi(t)}. 


4 . Numerical studies 

We conduct a number of simulation experiments to assess the performance of the pro¬ 
posed SN-based tests in comparison with the alternative methods in the literature. We 
generate functional data on a grid of 10 3 equispaced points in [ 0 , 1 ], and then convert dis¬ 
crete observations into functional objects by using B-splines with 20 basis functions. We 
also tried 40 and 100 basis functions and found that the number of basis functions does 
not affect our results much. Throughout the simulations, we set the number of Monte 
Carlo replications to be 1000 except for the i.i.d. bootstrap method in Benko et al. [1], 
where the number of replications is only 250 because of high computational cost. 

4.1. Comparison of covariance operators 

To investigate the finite sample properties of N (d) for dependent functional data, we 
modify the simulation setting considered in Panaretos et al. [16]. Formally, we consider 
the model, 


3 

'^2{Cj,iV2sm(2njt) + C j ,2S^cos{2njt)}, i = 1,2,.. ,,t G [0,1], (4.1) 

3 =1 

where the coefficients & = (£i 1,^2 1^3 1j£i,2>£2 2? £3 2/ are generated from a VAR pro¬ 
cess, 

& = + t/l - P 2 e;, (4.2) 

with e* € M 6 being a sequence of i.i.d. normal random variables with mean zero and co- 
variance matrix £ e = 1+ W diag(v) + We generate two independent functional 
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time series {Xi(t)} and {Y)(t)} from (4.1) with p = 0.5 and p = l. We compare the SN- 
based test with the PKM test which is designed for independent Gaussian process, and 
the traditional test which is constructed based on a consistent LRV estimator (denoted 
by CLRV), that is, GcL,N{d) = NctN^^^N, where S a is a lag window LRV estimator 
with Bartlett kernel and data dependent bandwidth (see Andrews (1991)). We report the 
simulation results for A T i = W = 100,200, K = 1,2,3,4,5 (c? = 1,3,6,10,15) and various 
values of v in Table 1. Results in scenario A show that the size distortion of all the three 
tests increases as K gets larger. The SN-based test has the best size compared to the 
other two tests. The PKM test is severely oversized due to the fact that it does not take 
the dependence into account. It is seen from the table that the CLRV test also has severe 


Table 1. Empirical sizes and size-adjusted powers of (i) the SN-based test, (ii) the PKM test 
and (iii) the CLRV test for testing the equality of the covariance operators. The nominal level 
is 5% 



Parameter 

W =n 2 


K 

1 

2 

3 

4 

5 

Kt 

Ki 

A 

v x 

= (12,7,0.5,9,5,0.3) 

100 

(i) 

4.3 

5.7 

6.8 

8.7 

14.3 

8.7 

10.7 





(ii) 

14.5 

20.9 

22.9 

32.2 

39.5 

32.0 

34.0 





(iii) 

9.1 

12.9 

20.9 

39.8 

67.9 

38.8 

48.6 


v Y 

= (12,7,0.5,9,5,0.3) 

200 

(i) 

4.7 

5.7 

4.6 

7.0 

8.0 

7.0 

7.0 





(ii) 

12.8 

20.6 

26.7 

34.7 

42.6 

34.8 

37.3 





(iii) 

6.9 

9.6 

14.5 

25.2 

41.5 

25.1 

28.9 

B 

v x 

= (14,7,0.5,6,5,0.3) 

100 

(i) 

19.1 

23.6 

17.7 

14.2 

12.7 

14.1 

13.0 





(ii) 

27.6 

37.7 

31.6 

22.9 

21.2 

23.1 

22.7 





(iii) 

27.0 

33.8 

23.0 

20.5 

14.0 

20.5 

15.8 


v Y 

= (8,7,0.5,6,5,0.3) 

200 

(i) 

31.2 

37.7 

30.4 

21.9 

21.9 

21.9 

22.1 





(ii) 

39.1 

61.6 

51.7 

44.2 

41.2 

44.2 

41.5 





(iii) 

37.6 

57.0 

44.7 

30.1 

24.3 

30.1 

24.9 

C 

v x 

= (12,7,0.5,9,3,0.3) 

100 

(i) 

5.5 

10.9 

30.8 

62.4 

64.7 

62.3 

63.7 





(ii) 

4.7 

16.1 

57.3 

94.6 

98.7 

94.4 

97.1 





(iii) 

5.5 

13.4 

42.3 

79.4 

70.8 

79.2 

74.3 


v Y 

= (12,7,0.5,3,9,0.3) 

200 

(i) 

5.3 

10.0 

45.7 

90.3 

94.3 

90.4 

92.5 





(ii) 

6.4 

13.0 

67.8 

99.9 

100.0 

99.9 

100.0 





(iii) 

6.1 

12.5 

60.5 

99.9 

99.8 

99.9 

99.8 

D 

v x 

= (12,7,0.5,9,5,0.3) 

100 

(i) 

6.1 

8.3 

28.3 

80.1 

82.2 

75.6 

80.6 





(ii) 

5.5 

14.6 

47.2 

100.0 

100.0 

94.7 

100.0 





(iii) 

6.9 

12.3 

37.2 

95.7 

88.6 

90.6 

90.7 


vy 

= (12,7,0.5,0,5,0.3) 

200 

(i) 

5.7 

8.9 

39.7 

96.3 

98.4 

95.5 

98.3 





(ii) 

6.4 

14.5 

53.6 

100.0 

100.0 

99.4 

100.0 





(iii) 

6.0 

12.9 

47.7 

100.0 

100.0 

99.3 

100.0 


Note: Under the alternatives, we simulate the size-adjusted critical values by assuming that both {Xi} 
and {Yi} are generated from (4.1) with p = 0.5, // = I and v = vx ■ 
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size distortion especially for large K , which is presumably due to the poor estimation 
of the LRV matrix of dx when the dimension is high. Under the alternatives, we report 
the size-adjusted power which is computed using finite sample critical values based on 
the simulation under the null model where we assume that both {W(i)} and {Yj(i)} are 
generated from (4.1) with p = 0.5, p = 1 and v = vj. From scenarios B-D in Table 1, 
we observe that the PKM is most powerful which is largely due to its severe upward size 
distortion. The SN-based test is less powerful compared to the other two tests but the 
power loss is generally moderate in most cases. Furthermore, we present the results when 
choosing K by 


K* = argminj 1 < J < 20: >a*\, j = 1,2, (4.3) 

L £i=i A xy J 

where ay = 85% and a 2 = 95%. An alternative way of choosing AT is to consider the penal¬ 
ized fit criteria (see Panaretos et al. [16] for the details). We notice that the performance 
of all the three tests based on automatic choice AT* is fairly close to the performance when 
AT = 4 or 5 in most cases. To sum up, the SN-based test provides the best size under 
the null and has reasonable power under different alternatives considered here, which is 
consistent with the “better size but less power” phenomenon seen in the univariate setup 
(Lobato [14] and Shao [21]). 

4.2. Comparison of eigenvalues and eigenfunctions 

In this subsection, we study the finite sample performance of the SN-based test for test¬ 
ing the equality of the eigenvalues and eigenfunctions. We consider the data generating 
process, 

2 

^2{Cj,iV2sm(2njt + Sj) + Cj, 2 V2cos(2njt + Sj)}, i = l,2,...,te [0,1], (4.4) 

j =i 

where £* = 21 ^ 2 , 2 )' i s a 4-variate VAR process (4.2) with e t £ R 4 being 

a sequence of i.i.d. normal random variables with mean zero and covariance matrix 
E e = diag(v) + I 4 I 4 . We set p = 0.5 and p = 0. Under H 2 , 0 (or if 2 ,a), {-^(t)} 

and {Yi(t)} are generated independently from (4.4) with 5\ = S 2 = 0 and =vy (or 
v x 7 Wy). Notice that the eigenvalues of (Aj(t)} and {V(t)} are given respectively, by 
vx and vy when <5i = S 2 = 0. Under Ff 3i0 and H 3 ^, we generate {Xi(t)} and {Y)(t)} 
independently from (4.4) with vx = vy, 5x, 1 — 5y, 1 = 5 , and Ay ,2 = Ay ,2 = 0, where <5 = 0 
under the null and 6 ^ 0 under the alternatives. We aim to test the equality of the first 
four eigenvalues and eigenfunctions separately and jointly. Because functional data are 
finite dimensional, we implement the untrimmed version of the SN-based tests, that is, 
e = 0. To further assess the performance of the SN-based test, we compare our method 
with the subsampling approach with several choices of subsampling widths and the i.i.d. 
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bootstrap method in Benko et al. [1]. Suppose N\ = N 2 = Ng. Let l be the subsam¬ 


pling width and X J suh i = X J S . 


sub.X.i 


- A) 


sub.y.i’ 


i = 1,2,...,sjv o (0 = Wl, where X J S 


sub.X.i 


and Ag ub Yi are estimates of the j th eigenvalues based on the ith nonoverlapping sub¬ 
samples and {Yk(t)}^ =( ^ i _ 1 ^ l+1 , respectively. The subsampling variance 

estimate is given by a 2 uhj = Y,Zf\Kub,i ~ Ei=° ( ° febyfe and the test 

statistic based on the subsampling variance estimate for testing the equality of the jth 
eigenvalue is defined as G su b,jv = Nq(X j x Nq — X Y n 0 ) 2 I a snh j- Since the data-dependent 
rule for choosing the subsampling width is not available in the current setting, we tried 
l = 8,12,16 for Nq = 48,96. For testing the equality of eigenvalues jointly and equality 
of the eigenfunctions, we shall consider a multivariate version of the subsampling-based 
test statistic which can be defined in a similar fashion. Table 2 summarizes some selective 
simulation results for testing the eigenvalues with various values of v. From scenario A, 
we see that performance of the SN-based test under the null is satisfactory while the size 
distortion of the subsampling-based method is quite severe and is sensitive to the choice 
of block size l. It is also not surprising to see that the i.i.d. bootstrap method has obvious 
size distortion as it does not take the dependence into account. Under the alternatives 
(scenarios B D), we report the size-adjusted power by using the simulated critical values 
as described in previous subsection. When the sample size is 48, the SN-based method 
delivers the highest power among the tests and it tends to have some moderate power 
loss when the sample size increases to 96. On the other hand, the subsampling method 
is sensitive to the choice of subsampling width and its power tends to decrease when a 
larger subsampling width is chosen. 

To test the equality of the first four eigenfunctions, we implement the SN-based test 
and the subsampling-based test with the basis functions, 


= (fey + feyj ■ • ■, fey + fe 


fey + fey > • • ■) fey 


l 'XY > 


v XY)i 


1 <j< 4,p = 4, 


(4.5) 


for testing individual eigenfunction and 


>r = (fey 


fey + fey • ■ • > fey 


fey)> 


1 < j <j*,p = 4, 


(4.6) 


with j* = 2,3,4, for testing the first j* eigenfunctions jointly (correspondingly Mg = 
3,5,6). The tests with the above basis functions tend to provide similar sizes but higher 
powers as compared to the tests with the basis functions zA in our simulation study. The 
basis functions v* is constructed by adding the same estimated eigenfunction (f> J XY to 
each component of , and the associated SN-based test is expected to be asymptoti¬ 
cally valid in view of the Bahadur representation (2.9). Selective simulation results are 
summarized in Table 3 and Figure 1 which present the sizes of the SN-based test, the 
subsampling-based test and the i.i.d. bootstrap method, and the size adjusted powers of 
the former two respectively. It is seen from Table 3 that the sizes of the SN-based test 
are accurate while the subsampling-based test is apparently size-distorted. It is some¬ 
what surprising to see that the i.i.d. bootstrap provides better sizes compared to the 
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Table 2. Empirical sizes and size-adjusted powers of (i) the SN-based test, the subsampling- 
based test with (ii) l = 8, (iii) l = 12 and (iv) l = 16, and (v) Benko et al.’s i.i.d. bootstrap 
based method for testing the equality of the first two eigenvalues separately (the columns with 
M = 1,2) and jointly (the column with M = (1,2)), and the equality of the first four eigenvalues 
jointly (the column with M = (1,2,3,4)). The nominal level is 5% and the number of replications 
for i.i.d. bootstrap method is 250 



Parameter 

£ 

II 

£ 


M 

1 

2 

(1, 2) 

(1, 2, 3, 4) 

A 

vx = (10,0.5,5,0.3) 

48 

(i) 

5.4 

5.1 

4.6 

3.8 




(ii) 

24.2 

38.5 

52.4 

90.8 




(iii) 

21.9 

28.8 

51.3 

68.8 




(iv) 

21.8 

28.1 

57.9 

44.7 




(v) 

11.2 

9.2 

11.6 

11.6 


vy = (10,0.5,5,0.3) 

96 

(i) 

5.2 

5.6 

4.8 

5.1 




(ii) 

19.0 

40.4 

46.4 

84.8 




(iii) 

16.3 

29.6 

38.2 

77.0 




(iv) 

16.0 

25.3 

36.5 

78.5 




(v) 

14.4 

8.4 

15.2 

15.2 

B 

vx = (20,0.5,5,0.3) 

48 

(i) 

25.1 

4.3 

21.8 

15.5 




(ii) 

24.2 

5.4 

13.3 

7.1 




(iii) 

19.8 

6.8 

00 

oo 

00 

oo 




(iv) 

14.1 

6.8 

8.0 

9.1 


vy = (10,0.5,5,0.3) 

96 

(i) 

48.4 

4.8 

35.6 

25.0 




(ii) 

58.4 

6.9 

29.4 

11.4 




(iii) 

50.9 

6.1 

29.7 

13.3 




(iv) 

53.8 

6.0 

29.3 

11.4 

C 

vx = (10,0.5,5,0.3) 

48 

(i) 

6.2 

70.6 

58.9 

44.1 




(ii) 

5.5 

68.1 

54.6 

13.9 




(iii) 

4.8 

49.3 

23.0 

16.9 




(iv) 

6.1 

34.2 

15.4 

21.2 


vy = (10,0.5,1,0.3) 

96 

(i) 

4.7 

91.4 

84.6 

77.6 




(ii) 

4.7 

98.7 

96.3 

69.7 




(iii) 

5.5 

97.9 

92.5 

51.6 




(iv) 

5.4 

96.5 

83.0 

29.4 

D 

v x = (20,0.5,5,0.3) 

48 

(i) 

27.0 

70.1 

68.4 

55.3 




(ii) 

25.8 

65.7 

51.1 

14.0 




(iii) 

20.9 

58.4 

23.9 

19.9 




(iv) 

14.9 

40.1 

11.8 

17.2 


vy = (10,0.5,1,0.3) 

96 

(i) 

55.3 

87.9 

88.4 

83.3 




(ii) 

54.5 

98.3 

96.9 

62.3 




(iii) 

48.8 

97.6 

95.2 

53.2 




(iv) 

50.1 

95.2 

88.0 

27.7 


Note: Under the alternatives, we simulate the size-adjusted critical values by assuming that both { X ,} 
and {11} are generated from (4.1) with p = 0.5, p, = 0 and v = vx ■ 
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Table 3. Empirical sizes of (i) the SN-based test, the subsampling-based test with (ii) 1 = 8, 
(iii) l = 12 and (iv) l = 16, and (v) Benko et al.’s i.i.d. bootstrap based method for testing the 
equality of the first two eigenfunctions separately (the columns with M = 1,2) and jointly (the 
column with M = (1,2)), and the equality of the first four eigenfunctions jointly (the column 


with M = (1,2,3,4)). The nominal level is 
is 250 

5% and the number of replicat 

ions for 

i.i.d. bootstrap 

Parameter 

£ 

II 

£ 


M 

1 

2 

(1. 2) 

(1, 2, 3, 4) 

A v x = (10,0.5,5,0.3) 

48 

(i) 

6.4 

3.2 

4.9 

4.8 



(h) 

39.5 

36.0 

78.7 

67.0 



(hi) 

62.9 

62.3 

24.8 

26.4 



(iv) 

32.6 

25.9 

9.6 

7.7 



(v) 

2.4 

11.2 

2.4 

2.4 

■vy = (8,0.5,4,0.3) 

96 

(i) 

4.5 

3.3 

4.3 

4.7 



(ii) 

18.2 

16.8 

27.2 

45.8 



(hi) 

24.9 

21.0 

49.0 

71.2 



(iv) 

32.6 

30.6 

75.9 

61.0 



(v) 

3.2 

12.4 

4.8 

6.8 

B yx = (4,0.5,2,0.3) 

48 

(i) 

8.2 

3.8 

7.0 

6.1 



(ii) 

43.3 

45.8 

83.4 

71.4 



(hi) 

66.6 

65.2 

23.8 

18.1 



(iv) 

26.5 

22.5 

5.8 

3.7 



(v) 

2.4 

6.0 

3.6 

1.6 

vy — (2,0.5,1,0.3) 

96 

(i) 

5.3 

4.5 

5.2 

4.9 



(ii) 

20.3 

24.2 

36.7 

54.6 



(hi) 

25.3 

27.9 

53.0 

75.7 



(iv) 

33.1 

32.6 

78.1 

60.1 



(v) 

2.8 

8.4 

6.8 

3.6 


subsampling-based approach which is designed for dependent data. Figure 1 plots the 
(size-adjusted) power functions of the SN-based test and the subsampling-based test 
which are monotonically increasing on <5. When N± = N 2 = 48, the SN-based test delivers 
the highest power in most cases. The subsampling-based test with a small subsampling 
width becomes most powerful when sample size increases to 96. Overall, the SN-based 
test is preferable as it provides quite accurate size under the null and has respectable 
power under the alternatives. 


5. Climate projections analysis 

We apply the SN-based test to a gridded spatio-temporal temperature dataset covering a 
subregion of North America. The dataset comes from two separate sources: gridded ob¬ 
servations generated from interpolation of station records (HadCRU), and gridded simu- 
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M=(1,2,3,4); N=48 



M=(1,2,3,4); N=96 



delta 


delta 


Figure 1. Size-adjusted powers of the SN-based test and the subsampling-based tests for testing 
the equality of the first two eigenfunctions separately and jointly, and the equality of the first 
four eigenfunctions jointly. 
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P-values 



120°W 110°W 100°W 90°W 80“W 


0 1 2 3 4 5 

Figure 2. p-values for testing the nullity of lag zero cross-correlation between the station obser¬ 
vations and model outputs at each location. The numbers 0-5 denote the ranges of the p-values, 
that is, 0 denotes [0.1, 1]; 1 denotes [0.05, 0.1]; 2 denotes [0.025, 0.05]; 3 denotes [0.01, 0.025]; 4 
denotes [0.005, 0.01] and 5 denotes [0, 0.005]. 


lations generated by an AOGCM (NOAA GFDL CM2.1). Both datasets provide monthly 
average temperature for the same 19-year period, 1980-1998. Each surface is viewed as a 
two-dimensional functional datum. The yearly average data have been recently analyzed 
in Zhang et al. [23], where the goal is to detect a possible change point of the bias be¬ 
tween the station observations and model outputs. In this paper, we analyze the monthly 
data from 1980 to 1998, which includes 228 functional images for each sequence. We fo¬ 
cus on the second-order properties and aim to test the equality of the eigenvalues and 
eigenfunctions of the station observations and model outputs. To perform the analysis, 
we first remove the seasonal mean functions from the two functional sequences. At each 
location, we have two time series from the demeaned functional sequences. We apply the 
SN-based test developed in Shao [21] to test whether their cross-correlation at lag zero 
is equal to zero. The p-values of these tests are plotted in Figure 2. The result tends to 
suggest that the dependence between the station observations and model outputs may 
not be negligible at certain regions as the corresponding p-values are extremely small. 
The two sample tests introduced in this paper are useful in this case because they are 
robust to such dependence. 

We perform FPCA on the demeaned sequences. Figure 3 plots the first three PC’s of the 
station observations and model outputs. We then apply the SN-based tests Ggjy N (M) 

and Gg 3 ^- N {M 0 ) (with p = 3) to the demeaned sequences, which yields the results sum¬ 
marized in Table 4. It is seen from the table that the first two eigenvalues of the station 
observations and model outputs may be the same, at least statistical significance is below 
the 10% level, while there is a significant difference between their third eigenvalue. The 
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PC3 (0.306, 9.5%) 



h - i r > 


120°W 110°W 100°W 90°W 80°W 



-2°C -1°C 0°C 1°C 2°C 


PC3 (0.548, 13.3%) 


1 — 



120°W 110°W 100°W 90°W 80°W 



-2°C -1°C 0°C 1°C 2°C 


Figure 3. The first three PCs of the station observations (left panels) and model outputs (right 
panels), and the associated eigenvalues and percentage of variations explained. 
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Table 4. Comparison of the eigenvalues and eigenfunctions of the covariance operators of the 
station observations and model outputs 


K 

Gsn,n(M) 

p -value 

o) 

p -value 

1 

10.8 

(0.1,1) 

126.4 

(0.025,0.05) 

2 

5.4 

(0.1,1) 

295.4 

(0,0.005) 

3 

119.9 

(0.005,0.01) 

34.2 

(0.1,1) 

- 

326.2 

(0.005,0.01) 

318.0 

(0.005,0.01) 


Note: The first three rows show the results for testing individual eigencomponent, and the last row shows 
the results for testing the first three eigenvalues and eigenfunctions jointly. 


SN-based tests also suggest that there are significant differences of the first and second 
PCs between the station observations and model outputs as the corresponding p -values 
are less than 5% while the difference between the third PCs is not significant at the 10% 

level; compare Figure 3. We also tried the basis functions v* and u** for G y s ^ N (M 0 ) 
(see (4.5) and (4.6)), which leads to the same conclusion. To sum up, our results suggest 
that the second-order properties of the station observations and model outputs may not 
be the same. 

In climate projection studies, the use of numerical models outputs has become quite 
common nowadays because of advances in computing power and efficient numerical al¬ 
gorithms. As mentioned in Jun et al. [11], “Climate models are evaluated on how well 
they simulate the current mean climate state, how they can reproduce the observed cli¬ 
mate change over the last century, how well they simulate specific processes, and how 
well they agree with proxy data for very different time periods in the past.” Further¬ 
more, different institutions produce different model outputs based on different choices of 
parametrizations, model components, as well as initial and boundary conditions. Thus 
there is a critical need to assess the discrepancy/similarity between numerical model 
outputs and real observations, as well as among various model outputs. The two sample 
tests proposed here can be used towards this assessment at a preliminary stage to get 
a quantitative idea of the difference, followed by a detailed statistical characterization 
using sophisticated spatio-temporal modeling techniques (see, e.g., Jun et al. [11]). In 
particular, the observed significance level for each test can be used as a similarity index 
that measures the similarity between numerical model outputs and real observations, and 
may be used to rank model outputs. A detailed study along this line would be interesting, 
but is beyond the scope of this article. 
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Supplement to “Two sample inference for the second-order property of tem¬ 
porally dependent functional data” (DOI: 10.3150/13-BEJ592SUPP; .pdf). This 
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